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ABSTRACT 

We propose easy ways of correcting for the systematic errors caused by the photon 
noise and the pixelation effect in cosmic shear measurements. Our treatment of noise 
can reliably remove the noise contamination to the cosmic shear even when the flux 
density of the noise is comparable with those of the sources. For pixelated images, 
we find that one can accurately reconstruct their corresponding continuous images 
by interpolating the logarithms of the pixel readouts with either the Bicubic or the 
Bicubic Spline method as long as the pixel size is about less than the scale size of the 
point spread function (PSF, including the pixel response function), a condition which 
is almost always satisfied in practice. Our methodology is well defined regardless of 
the morphologies of the galaxies and the PSF. Despite that our discussion is based 
on the shear measurement method of Zhang (2008), our way of treating the noise can 
in principle be considered in other methods, and the interpolation method that we 
introduce for reconstructing continuous images from pixelated ones is generally useful 
for digital image processing of all purposes. 

Key words: cosmology: gravitational lensing - methods: data analysis - techniques: 
image processing: large scale structure 



1 INTRODUCTION 

Weak gravitational lensing has been widely used as a direct 
probe of the large scale structure (see reviews by Bartelmann 
& Schneider 2001; Wittman 2002; Refregier 2003). By mea- 
suring the systematic distortions of background galaxy im- 
ages, one can place constraints on the cosmological param- 
eters (Bacon et al. 2000; Kaiser et al. 2000; van Waerbeke 
et al. 2000; Wittman et al. 2000; Maoli et al. 2001; Rhodes 
et al. 2001; van Waerbeke et al. 2001; Hoekstra et al. 2002; 
Refregier et al. 2002; Bacon et al. 2003; Brown et al. 2003; 
Hamana et al. 2003; Jarvis et al. 2003; Rhodes et al. 2004; 
Heymans et al. 2005; Massey et al. 2005; van Waerbeke et 
al. 2005; Dahle 2006; Hoekstra et al. 2006; Jarvis et al. 2006; 
Semboloni et al. 2006; Hetterscheidt et al. 2007; Schrabback 
et al. 2007). With accurate redshift information, the geom- 
etry and the structure growth rate of our Universe can be 
constrained as functions of redshift separately, providing a 
consistency test of the gravity theory (Hu 2002; Abazajian & 
Dodelson 2003; Jain & Taylor 2003; Acquaviva et al. 2004; 
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Bernstein & Jain 2004; Hu & Jain 2004; Kratochvil et al. 
2004; Song & Knox 2004; Takada & Jain 2004; Takada & 
White 2004; Ishak 2005; Simpson & Bridle 2005; Song 2005; 
Zhang et al. 2005; Hannestad et al. 2006; Ishak et al. 2006; 
Zhan 2006; Knox et al. 2006; Schimd et al. 2007; Taylor et 
al. 2007). 

A key issue in weak lensing is about how to measure 
the cosmic shear with galaxy shapes. This is difficult mainly 
because the signal-to-noise ratio of the measurement on one 
galaxy is typically only a few percent. It is therefore ex- 
tremely important for any shear measurement method to 
carefully treat any possible systematic errors, including at 
least the following: the correction due to the image smearing 
by the point spread function (PSF, including the pixel re- 
sponse function); the photon noise; the pixelation effect due 
to the discrete nature of the CCD pixels. There have been 
a number of methods proposed to deal with the corrections 
due to the PSF (see Tyson et al. 1990; Bonnet & Mellier 
1995; Kaiser et al. 1995; Luppino & Kaiser 1997; Hoekstra et 
al. 1998; Rhodes et al. 2000; Kaiser 2000; Bridle et al. 2001; 
Bernstein & Jarvis 2002; Refregier & Bacon 2003; Massey & 
Refregier 2005; Kuijken 2006; Miller et al. 2007; Nakajima 
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& Bernstein 2007; Kitching et al. 2008; Zhang 2008). How- 
ever, the photon noise and the pixelation effect remain to 
be treated in a more systematic way. For example, existing 
model fitting methods {e.g. , Bridle et al. 2001; Bernstein 
& Jarvis 2002; Kuijken 2006; Miller et al. 2007; Nakajima 
& Bernstein 2007; Kitching et al. 2008) use the variance of 
the noise to weight the pixels in their chi-square fittings. 
It is not clear to what level the noise contamination to the 
shear recovery can be removed in this way, especially for 
correlated background noise (see, e.g. , Massey et al. 2007). 
By integrating the model over the pixels, the model fitting 
methods essentially use linear interpolations to treat the pix- 
elation effect. This is found not accurate as we will discuss 
in §3.2. Indeed, as we will discuss later in the paper, there 
are two types of noise: the astronomical photon noise and 
the "photon counting" shot noise. While the second type 
diminishes when the exposure time increases, the first type 
does not. We will mainly focus on the astronomical noise 
in this paper. For the counting shot noise, we will simply 
argue in §5 that its contamination to the shear recovery can 
be significantly suppressed by increasing the exposure time. 
Since these issues are not specifically addressed in any pre- 
vious weak lensing literatures, it is important to point them 
out in this paper. 

In a recent work by Zhang (2008) (Z08 hereafter) , a new 
and simple way of measuring the cosmic shear is found. Its 
main advantages includes: 1. it is mathematically simple; 2. 
it is free of assumptions on the morphologies of the galaxies 
and the PSF; 3. it enables us to probe the shear information 
from galaxy substructures, thereby improving the signal-to- 
noise ratio. These facts encourage us to extend the method 
further by including a treatment of the photon noise and 
the pixelation effect. Fortunately, we find that these two 
types of systematic errors can be treated in a simple and 
model-independent way based on the method of Z08. As 
will become clear later in this paper, the method we adopt 
to remove the noise contamination can also be considered 
in other shear measurement methods, and our treatment of 
the pixelation effect is generally useful for image processing 
of all purposes. 

The paper is organized as follows: in §2, we briefly re- 
view the shear measurement method of Z08; in §3, we show 
how to treat the photon noise (in §3.1) and the pixelation 
effect (in §3.2) in weak lensing; in §4, we use computer- 
generated mock galaxy images to test the performance of 
our method; finally, we summarize and discuss remaining 
issues in §5. 



where A^- = Sij + <&ij, and <E>;j = dS9{ /d6f are the spatial 
derivatives of the lensing deflection angle. Matrix A is often 
expressed in terms of the convergence k — (& xx + & yy )/2 
and the two shear components 71 = (& xx — & yy )/2 and 
72 = &x y . Assuming the intrinsic galaxy image fs(6 S ) is 
statistically isotropic, the shear components can be simply 
related to the derivatives of the surface brightness field as 
(Seljak & Zaldarriaga 1999): 



1 {{d*fi 
7i = -7, 



(dyflf) 



72 = — 



2 ((ft,//) 2 + (c^) 2 ) 
(d x fidyfi) 

{(d x fjy + (c\/,) 2 > 



(2) 



where the averages are taken over the galaxy. 

Eq.[2] is useful only when the angular resolution of 
the observation is infinitely high. In practice, the observed 
galaxy surface brightness distribution fo is always equal to 
the lensed galaxy image // convoluted with the PSF, i.e. : 



fo(0) 



(3) 



where W is the PSF. Z08 has shown how to modify eq.(2) 
when the PSF is an isotropic Gaussian function, which can 
be written as: 



(4) 



where j3 is the scale radius of the Gaussian function. The new 
relation between the shear components and the derivatives 
of the surface brightness field is: 



71 = 

72 = 
where 



1 ({d x fof - (dyfp) 2 ) 
2{(d x fo) 2 + (d y fo) 2 ) + A 

(d x fod y fo) 
~{(d x fo) 2 + (d y fo) 2 )+A 



-(V/o-V(V 2 /o)) 



(5) 



(6) 



For a general PSF, one can transform it into the desired 
isotropic Gaussian form through a convolution in Fourier 
space. The scale radius /3 of the target PSF should be 
larger than that of the original PSF to avoid singularities in 
the convolution. Furthermore, as shown in Z08, the spatial 
derivatives required by eq.(5) can also be easily evaluated in 
Fourier space. 



2 REVIEW OF THE SHEAR MEASUREMENT 
METHOD 

Z08 proposes a way of measuring the cosmic shear with the 
spatial derivatives of the galaxy surface brightness field. To 
do so, let us define the surface brightness on the image plane 
as fiiO 1 ), and that on the source plane as fs{0 s ), where 9 1 
and 9 s are the position angles on the image and source plane 
respectively. These quantities are related in a simple way as: 

fifi 1 ) = fs(0 S ) (1) 



3 TREATING THE PHOTON NOISE AND 
THE PIXELATION EFFECT 

In this section, we introduce the basic ideas for treating the 
photon noise and the pixelation effect in §3.1 and §3.2 re- 
spectively. Numerical examples are given in the next section. 

3.1 Photon Noise 

First of all, there are two types of photon noise: the astro- 
nomical photon noise due to the fluctuation of the back- 
ground, and the "photon counting" shot noise due to finite 
exposure time. Note that the first type of noise, like the 
source, is convoluted by the PSF, while the second type 
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of noise varies from pixel to pixel even if the pixel size is 
much smaller than the PSF size. In the rest of the paper, 
we mainly deal with the astronomical noise. For the "pho- 
ton counting" shot noise, we will simply argue in §5 that 
its contamination to the shear recovery can be significantly 
suppressed by increasing the exposure time. 

The presence of the photon noise makes the measure- 
ment of the cosmic shear more complicated in two ways: 1. 
the observed surface brightness fo is from both the lensed 
source and the un-lensed foreground noise; 2. because of the 
aliasing power caused by the non-periodic boundaries of the 
noisy map, the measurement of the spatial derivatives of the 
surface brightness field cannot be accurately performed in 
Fourier space. Note that simple treatments such as filtering 
out the noise outside of the source image do not completely 
fix this problem, because the noise inside the image can still 
bias the shear estimate. Fortunately, as we show in the rest 
of this section, our master equation [eq.(5)] for estimating 
the cosmic shear can be easily adapted to solve both prob- 
lems. 

In the method of Z08, to isolate the source signals in a 
noisy map, let us first write the total observed surface bright- 
ness fo as the sum of the contributions from the source f 
and the noise f n , i.e. , fo = f + f n - Note that in this 
case, instead of fo, f s should be used in eq.(5) to correctly 
measure the shear components. For this purpose, let us use 
the following relation: 

{(d*ff) = ((d x fo) 2 )-{(d x .n 2 ) (7) 

- 2(d x fd x f n ) 

{(dyf) 2 ) = {(dyfof) - {(dyf"f) 

- 2(d y fd y n 
(d.rd y n = (d*fod y f )-(d x f n d y n 

- (d y f s d x n - (d x fd y n 
(vr-v(v 2 r)> = <v/o • v(v 2 /o)) - <vr • v(v 2 d) 

- (vr •v(v 2 r)>-<vr-v(v 2 r)> 

For simplicity, in this paper, we only consider the photon 
noise that is from the foreground or the instruments. The 
surface brightness distribution of the noise is therefore un- 
corrected with that of the background sources 1 . Under this 
assumption, the cross-correlations between the source and 
the noise terms (such as, e.g. , (d x f s d x f n )) should vanish. 
Eq.(7) therefore becomes: 

{(d x f) 2 ) = {(d x fo) 2 ) - {(d x f n f) (8) 

{(dyf) 2 ) = {(dyfof) - {(dyD 2 ) 

(dxf s d y f s ) = (d x fod y fo) - (d x f n d y f n ) 

<vr-v(v 2 r)> = <v/o • v(v 2 / )) - (vr • v(v 2 D) 

The relations in eq.(8) suggest an easy way of removing 
the contaminations from the photon noise: one can use a 
neighboring map of pure noise f n to estimate {{d x f n ) 2 }, 
({dyD 2 ), (9xf n d y f n ), (V/ n -V(V 2 / n )),and subtract them 
from their counterparts evaluated from the noisy source map 
fo to get the source terms required by eq.(5). Note that 
since the noise photons are distributed differently in each 

1 More general cases (e.g. , photon noise coming from faint back- 
ground sources) are more complicated, and will be dealt with in 
a future work. 



map, the above procedure does not exactly remove the noise 
contribution for each source image. However, the method is 
statistically accurate as long as the statistical properties of 
the photon noise are stable over a reasonably large scale. In 
other words, the differences in the noise distributions of two 
maps add statistical errors to the measured cosmic shear 
through this procedure, but no systematic errors. Finally, 
the pure noise map should be a close neighbor of the source 
map so that they share the same point spread function. 

To evaluate the derivatives of the surface brightness 
field of a noisy map in the Fourier space, we need to deal with 
the non-periodic boundaries appropriately to avoid aliasing 
powers. This can be done by gradually attenuating the noise 
towards the boundaries of the map. The attenuation can 
take an arbitrary form as long as the following criterions are 
satisfied: 1. the source region is not affected; 2. the edges of 
the map should be rendered sufficiently faint; 3. the atten- 
uation amplitude should not have abrupt spatial variations; 
4. to properly remove the noise contamination, the same at- 
tenuation should also be applied to the neighboring map of 
pure noise. 

In §4, we show numerical examples to support the noise 
treatment discussed above. 

3.2 The Pixelation Effect 

Modern astronomical images are commonly recorded on 
CCD pixels, the discrete nature of which may affect the 
accuracy of the cosmic shear measurements. In general, to 
avoid significant shear measurement errors, the pixel size 
should be at least a few times smaller than the scale ra- 
dius (or FWHM) of the PSF. For instance, we find that the 
method of Z08 requires the scale radius of the PSF to be 
roughly 3 or 4 times larger than the pixel size. This require- 
ment is often not satisfied in space-based observations. It 
is therefore useful to have a method which can reconstruct 
continuous images from under-sampled ones. 

This is indeed a well-defined interpolation problem: how 
to reconstruct a continuous function if its value is only given 
at a set of discrete points. In the context of weak lensing, 
there is a quantitative way of testing the performance of 
the interpolation method, which is to check the accuracy of 
the shear recovery using the reconstructed images. The best 
method should yield the fastest convergence to an accurate 
image reconstruction (or shear recovery) as the pixel size 
becomes smaller. 

We have tested several standard 2D interpolation meth- 
ods, including Bilinear interpolation, Bicubic interpolation, 
and Bicubic Spline interpolation (see Press et al. 1992 for 
details) . Although these conventional methods perform rea- 
sonably well, we find that they can all be significantly im- 
proved by interpolating the natural logarithms of the data 
instead of the data themselves. This is mainly due to two 
reasons: 1. the values of the data have a lower bound — 
zero; 2. at large distances, the PSF typically falls off expo- 
nentially. For convenience, in the rest of the paper, we call 
such extensions of the three classical interpolation meth- 
ods their original names with the prefix "Log-" , and always 
abbreviate "Bicubic Spline" to "Spline". The mathemati- 
cal definitions of the six methods (i.e. , Bilinear, Bicubic, 
Spline, Log-Bilinear, Log-Bicubic, Log-Spline) are given in 
the appendix. In §4, we show that the Log-Bicubic and Log- 
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Spline methods are most accurate among the six. As will be 
shown in our numerical examples, continuous images that 
are reconstructed by these two interpolation methods yield 
negligible systematic errors in shear recovery as long as the 
pixel size is about smaller than the PSF size (twice its scale 
radius). Note that the pixel size is rarely larger than the 
PSF size in practice, because the pixel response function is 
a part of the PSF. 



4 NUMERICAL EXAMPLES 

We present numerical examples to support the ideas intro- 
duced in the previous section. The general setup of our 
numerical simulations are given in §4.1. In §4.2 and §4.3, 
we test our treatments of the photon noise and the pixela- 
tion effect separately. Finally, the overall performance of our 
method is shown in §4.4. 




Figure 1. Two PSF forms used in this paper. The letter on the 
up- left corner of each plot is the label of the PSF defined in eq.(9) 
(rotated by certain angles). The contours mark 0.0025%, 0.025%, 
0.25%, 2.5%, and 25% of the peak intensity. 



4.1 General Setup 

Each of our galaxy images is placed on a 2" x 2 n grid, where 
n is an integer (typical chosen to be 8). Note that such a 
choice facilitates the Fast Fourier Transformation (FFT). 
Each grid point is treated as the location of the center of a 
CCD pixel, whose side length is equal to the grid size. All 
sizes in our simulations are expressed in units of the grid 
size, i.e. , the pixel size. 

The form of the PSF is chosen from the following two 
functions rotated by certain angles: 



Wa{x, y) oc exp 
Wb{x, y) oc exp 
+ 0.03 exp 



2R 2 PSF 



2R PSF 



(x 2 + O.Sy 2 ) 
(x 2 + O.Sy 2 ) 



(9) 



E>2 



+ 0.2 



p2 



+ 0.2 



where Rpsf is the scale radius of the PSF. A schematic 
view of the PSF functions is shown in fig. 1 . Note that the 
additional term in Wb mimics the diffraction spikes. As dis- 
cussed in §2, before measuring the shear using eq.(5), the 
PSF is always transformed into the desired isotropic Gaus- 
sian form, whose scale radius /? should be slightly larger than 
Rpsf defined here to avoid singularities in the transforma- 
tion. Note that we have reserved the Greek letter /3 for the 
scale radius of the target PSF to distinguish it from Rpsf- 
Both the galaxies and the noise are treated as collections 
of point sources. For example, the image of a galaxy in our 
simulation is typically made of a few hundred or thousand 
points. The advantage of doing so is that one can easily lense 
the galaxy by displacing its point sources and modifying 
their amplitudes. The intensities of the point sources are dis- 
tributed to the neighboring grid points of their locations ac- 
cording to the PSF. For example, a point source of intensity 
A at location x contributes an intensity of A ■ WpsF{y — x) 
to the grid point at location y. The total intensity on a grid 
point is the sum of contributions from all the point sources. 
Since everything is composed of point sources in our sim- 
ulation, we will mostly call the surface brightness the flux 
density in the rest of the paper. 



There are two types of mock galaxies we use in this pa- 
per: regular disk galaxies and irregular galaxies. Our regular 
galaxy contains a thin circular disk of an exponential profile 
and a co-axial de Vaucouleurs-type bulge (de Vaucouleurs 
et al. 1991). On average, its face-on surface brightness dis- 
tribution is parameterized as: 

fir) oc exp(-r/r disk ) + f b/d exp [- [r/r bu i ge ) 1/4 ] 



(10) 



where r is the distance to the galaxy center, r bu i g e and rdisk 
are the scale radii of the bulge and the disk respectively, and 
f b / d determines the relative brightness of the bulge with re- 
spect to the disk. In the simulation, this profile is realized by 
properly and randomly placing a certain number (typically a 
few hundred) of point sources. These point sources are pro- 
jected onto a randomly oriented image plane, lensed, and 
finally assigned to the CCD pixels according to the PSF to 
yield the galaxy image. Our irregular galaxies are generated 
by 2D random walks. The random walk starts from the cen- 
ter of the grid, and continues for a certain number of steps. 
Each step has a fixed size and a completely random orien- 
tation in the image plane. The joint of every two adjacent 
steps gives the pre-lensing position of a point source of the 
galaxy in the image plane. The resulting irregular galaxies 
usually contain abundant substructures. 

For numerical manageability, we always cutoff the 
galaxy profile at a certain radius Rg, which is denoted as 
the scale radius of the galaxy. This is done by excluding the 
points that are outside of radius Ra in generating our reg- 
ular galaxies. For the irregular galaxies, the random walker 
is sent back to the origin to continue from there when it 
reaches the radius Ra- 

Without loss of generality, we set f b / d — 1/3, 
r d isk/r bu ige = 3, and r disk = Rg for the regular galaxies. 



4.2 Testing the Photon Noise Treatment 

As our first example, we use 10000 mock regular galaxies to 
test the treatment of photon noise discussed in §3.1. Each 
galaxy is made of 600 point sources. The galaxy size Rg is 
fixed at 30. The angle between the line-of-sight direction and 
the normal vector of the disk plane is randomly chosen from 
[0,tt/6]. The PSF we use is W A of eq.(9) with R PS f = 4 (no 
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pixelation effect). The scale radius ft of the target isotropic 
Gaussian PSF is 6. The noise is made of Poisson distributed 
point sources, the mean number density of which is roughly 
one per 4x4 area. Each point source of the noise is given the 
same intensity, the value of which can be freely adjusted to 
determine the ratio of the mean flux density of the noise to 
that of the galaxy (within its scale radius). Note that a sep- 
arate pure noise map is generated for each noisy galaxy map 
for removing the noise contaminations to the shear signal. 
Fig. 2 shows sample images of four different noise-to-galaxy 
flux density ratios: 0.02, 0.1, 0.5, 1.5. To avoid aliasing pow- 
ers in Fourier transformation, the noise near the boundaries 
of the map is filtered out by a window function defined as 
follows: 



WfiHer(r) 



1 

exp 



_ 1 ( r-fieore V 
2 ^ R a idth ) 



if T ^ Rcore: 
if T ^ Rcore ■ 



where r is the distance to the map/galaxy center, R core is 
the radius encircling the unaffected region, and Rmidth de- 
termines the width of the transition area. Note that the flat 
core of the filter should be sufficiently large to avoid affect- 
ing the galaxy images. In our simulations, we always choose 
Rcore = Rg + 4/?, and R wi dth = /3. Finally, to correct for 
the noise, for each galaxy map, we generate a map of pure 
noise to measure the noise properties required by eq.(8). The 
same filter is also applied to the pure noise map before the 
Fourier transformation. 

In fig.3, we plot the recovered shear values for different 
noise-to-galaxy flux density ratios In/Ig- The input cosmic 
shear (71,72) is (0.023,-0.037), displayed as dotted lines. 
The red data points with la error bars are our main re- 
sults achieved through the complete treatment of noise. For 
a comparison, the blue ones show the shear values measured 
directly from the noisy galaxy maps using eq.(5) without 
correcting for the noise (but the filter given by eq.(ll) is 
still applied to avoid the aliasing power in Fourier trans- 
formation). The figure clearly demonstrates that our noise 
treatment works remarkably well even when the mean flux 
density of the noise is comparable with that of the galaxy, 
though we caution that the size of the error bar grows with 
increasing noise intensity. On the other hand, the blue data 
points indicates that without a proper treatment of the 
noise, the measured shear values quickly drops to zero as the 
noise flux becomes dominant, deviating significantly from 
the input shear values. 

Finally, it is worth noting that both here and in the 
simulations reported in the rest of the paper, the number 
density of the noise points is always chosen to be roughly 
one point per PSF area (~ R% SF ) or slightly more. This is 
due to two reasons: 1. for a fixed mean noise flux density, 
a higher number density of the Poisson distributed noise 
points indeed leads to lower spatial fluctuations of the noise 
surface brightness field, therefore a less contamination to 
the shear signal, or a less challenging condition; 2. On the 
other hand, if the number density is much smaller than one 
point per PSF area, the image turns into a collection of 
discrete point-like sources, causing both the nominator and 
the denominator of eq.(5) to become small differences of 
large numbers, which are well known sources for numerical 
errors. The second point simply means it is hard to measure 
the shapes of sources that are much smaller than the size of 







Figure 2. Sample images of four different noise-to-galaxy flux 
density ratios. The up-left, up-right, lower-left, lower-right plots 
are for flux density ratios of 0.02, 0.1, 0.5, 1.5 respectively. 



the PSF. In practice, we can avoid the second situation by 
increasing the observation/integration time. 



4.3 Testing the Treatment of the Pixelation Effect 

When the scale radius of the PSF is less than 3 or 4 times the 
pixel size, the galaxy/noise images start to look pixelated, 
and the shear recovery accuracy may be strongly affected by 
the discrete nature of the CCD pixels. To reconstruct con- 
tinuous images, we use 2D interpolation methods to insert 
finer grid points. The finer grid size is chosen to be 2 m (m 
is an integer) times smaller than the original grid size, and 
at least less than a quarter of the PSF scale radius Rpsf- It 
is worth noting that interpolation of the pixelated image, if 
necessary, is always the first step in our shear measurement 
procedure. 

An example of a pixelated image is shown in the up- 
left corner of fig. 4, for which the scale radius of the PSF is 
0.5. The high resolution images reconstructed by the three 
interpolation methods and their logarithmic extensions dis- 
cussed in §3.2 are show in the lower half of fig.4. The true 
high resolution image is on the up-right corner of the fig- 
ure. By simply comparing the morphologies of the inter- 
polated images by eye, one may already tend to conclude 
that the performances of the three conventional methods 
are improved if we interpolate the log of the data instead 
of the data itself. For instance, negative intensities (denoted 
as blue regions in the figure) are commonly found in the 
interpolated maps by the Bicubic and Spline methods, but 
absent in those by the Log-Bicubic or Log-Spline methods; 
the filamentary features produced by the Bilinear method 
become somewhat less prominent in the map processed by 
the Log-Bilinear method. 
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Figure 3. The shear values with ltr error bars measured from im- 
ages of different noise-to-galaxy flux density ratios ijv/^G- The 
measurement uses 10000 mock regular galaxies. The dotted lines 
indicate the input shear values. The red data points are our main 
results achieved through the complete treatment of noise intro- 
duced in §3.1. The blue ones show the shear values measured di- 
rectly from the noisy galaxy maps using eq.(5) without removing 
the noise contaminations. 



Figure 5. The distribution of the shear values measured from 
the interpolated images of a single galaxy placed at random po- 
sitions on the grid. The purple, blue, red, and black histograms 
correspond to the pixel-size-to-PSF-radius ratio of 2, 1.5, 1, 0.5 
respectively. Each panel shows the results from a single interpo- 
lation method, whose name is indicated in the upper-left corner 
of the plot. All the histograms are normalized so that their peak 
values arc one. 
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Figure 4. The example of a pixelated image (upper-left), its high 
resolution counterpart (upper-right), and its interpolated images 
by the Bilinear, Bicubic, Spline, Log-Bilinear, Log-Bicubic, and 
Log-Spline methods (middle and lower panels). The blue regions 
have negative intensities. 



To test the interpolation methods more quantitatively, 
we may compare the shear values measured from the inter- 
polated maps using eq.(5) (in the absence of noise). Clearly, 
all the interpolation methods should yield the same and cor- 
rect shear estimates for a given galaxy and PSF if the PSF 
scale radius is much larger than the pixel size, i.e. , in the 
absence of the pixelation effect. For small PSF sizes, as we 
have seen, the continuous images interpolated by different 
methods look unlike each other, resulting in possibly very 
different shear values. Moreover, since the source is sparsely 
sampled in this case, the original (pixelated) image, the in- 
terpolated image, and the measured shear values all depend 
on the relative positions of the pixels with respect to the 
source. In fig. 5, we show the distributions (as histograms) 
of the shear values estimated from a single galaxy that is 
placed at different/random locations on the grid. For sim- 
plicity and clarity, we do not include any photon noise here. 
The ratio of the galaxy size Rg to the PSF scale radius Rpsf 
is fixed at 5. The ratio of the pixel size to Rpsf is chosen 
to be 2, 1.5, 1, 0.5, the results of which are represented by 
the purple, blue, red, and black histograms respectively. The 
figure shows that as the pixel size decreases relative to the 
PSF size, the shear distributions converge more rapidly to 
a delta function at the correct position in methods with the 
prefix "Log-" , manifesting again the value of the logarithmic 
extensions of the three classic interpolation methods. 

Finally, let us find out which interpolation method is 
best suited to weak lensing. For this purpose, we test the 
accuracy of shear recovery with a large number of inter- 
polated galaxy images. To make it a more convincing test, 
we use our morphologically rich irregular galaxies, each of 
which is generated by 1000 random steps. We consider three 
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choices for the random walk step size and the galaxy scale 
radius R G - {0.5R PS f, 6.67Rpsf), (0.25R PS f, 3.33R PS f), 
and (0.125-Rpsf, 1.67-Rpsf), referring to large, medium, 
and small galaxies respectively. The PSF we use is Wb ■ ft of 
the target isotropic Gaussian PSF is set to be 4/3 of Rpsf- 
The cosmic shear (71, 72) is chosen to be (—0.031, 0.018). 
No photon noise is included. Our results are summarized in 
fig. 6, in which we plot the measured shear values against 
the ratio of the pixel size to Rpsf- In the upper, middle, 
and lower panels of the figure, we report the results from 
averaging over 10000 large, medium, and small size galaxies 
respectively. The dotted lines refer to the input shear val- 
ues. The cyan, blue, magenta, green, red, and black data 
points with la error bars are from the Bilinear, Bicubic, 
Spline, Log-Bilinear, Log-Bicubic, and Log-Spline methods 
respectively. According to the figure, we can draw several 
conclusions: 

1. As a sanity check, we confirm that all the interpola- 
tion methods work well when the PSF size is much larger 
than the pixel size; 

2. Log-Bicubic and Log-Spline are the two most suc- 
cessful methods. Both of them can correctly recover the in- 
put shear as long as Rpsf is about larger than a half of 
the pixel size, regardless of the galaxy size. Note that one 
should not expect any interpolation method to work well 
when Rpsf <C 0.5 unless the source images are sufficiently 
smooth over the scale of the pixel size. 

3. The pixelation effect is less important for larger 
galaxies. For instance, by comparing the results in the three 
panels of fig. 6, we see that the quality of shear recovery 
becomes increasingly poor for galaxies of smaller sizes for 
a given Rpsf- This is not surprising because the struc- 
tures/shapes of large galaxies are better resolved than those 
of smaller ones. Meanwhile, it is encouraging to note that 
the Log-Bicubic and Log-Spline methods perform fairly well 
for Rpsf ij 0.5 even when the galaxy size is comparable to 
the PSF size. 



4.4 Testing the Overall Performance 

The purpose of this section is to test our shear measure- 
ment method under general conditions, i.e. , in the presence 
of both photon noise and the pixelation effect. Fig. 7 shows 
the pipeline of the numerical procedures we take in general 
cases. A detailed explanation of each item in the graph has 
been given in §3. 

In fig. 8, we show the shear recovery results for both reg- 
ular and irregular galaxies with different PSF forms. In each 
panel, we use 10000 mock galaxies and scale them to four 
different sizes (galaxy radius ranges from 2.5 to 10 times the 
scale radius of the PSF) to test the accuracy of shear recov- 
ery. In all panels, the PSF scale radius (Rpsf) is half of the 
pixel size, corresponding to roughly the maximum pixela- 
tion effect that can be treated by an interpolation method. 
The red and black data points are from the Log-Bicubic 
and Log-Spline methods respectively. The input shear val- 
ues are shown by dashed lines. The scale radius of the target 
isotropic Gaussian PSF is always 1.5Rpsf- To avoid alias- 
ing powers in Fourier transformation, we use eq.(ll) to filter 
the noise near the boundaries of the map. From the top to 
the bottom panel, the ratio of the mean flux density of the 
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Figure 6. The measured shear values plotted against the ratio of 
the pixel size to the PSF scale radius Rpsf- The upper, middle, 
and lower windows show the results from averaging over 10000 
relatively large, medium, and small size galaxies respectively. The 
definition of the galaxy size in this example can be found in §4.3. 
The cyan, blue, magenta, green, red, and black data points with 
la error bars are the results from the Bilinear, Bicubic, Spline, 
Log-Bilinear, Log-Bicubic, and Log-Spline interpolation methods 
respectively. The input shear values are shown as dotted lines. 



noise to that of the galaxy (In/Ig) is 0.1, 0.5, 0.6, and 0.2 
respectively. 

The figure indicates that our method generally works 
well on galaxies of sizes that are at least a few times larger 
than the PSF size. Small discrepancies between the input 
shear values and the measured ones do exist when the galaxy 
size is comparable to the PSF size. The residual systematic 
errors will be studied with a much larger ensemble of galaxies 
in another work. 
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Input the Galaxy Image and a 
Neighboring Map of Pure Noise 



Reconstruct Continuous (High Resolution) 
Images Using Interpolation 



Attenuate the Intensities on the Boundaries of 
Both Images with the Same Window Function 



Fourier Transform Both Images Using FFT 



Transform the PSF into an Isotropic 
Gaussian Form for Both Images 



Measure the Averages of the Spatial Derivatives of 
the Surface Brightness Fields of Both Images 



Correct for the Noise Contamination 



Calculate the Shear Components 
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Figure 7. The pipeline of our shear measurement procedures in 
the presence of both the photon noise and the pixelation effect. 

5 DISCUSSION AND CONCLUSIONS 

We have discussed how to correct for the systematic errors 
due to the photon noise and the pixelation effect in cosmic 
shear measurements. Our treatment of photon noise allows 
us to reliably remove the noise contamination to the cos- 
mic shear even when the noise flux density is comparable 
with that of the sources. In principle, our method works 
regardless of the brightness of the noise, though when the 
noise is much brighter than the sources, one needs to worry 
about image selections. To deal with pixelated images, our 
approach is to reconstruct continuous images by interpolat- 
ing the natural logarithms of the pixel readouts with either 
the Bicubic or Bicubic Spline method. This technique is ac- 
curate for the purpose of shear recovery as long as the scale 
radius of the PSF is larger than about a half of the pixel 
size, a condition which is almost always satisfied in practice. 

Despite the fact that our study has been based on the 
shear measurement method of Z08, a part of our methodol- 
ogy is generally useful for other shear measurement methods, 
or even other astronomical measurements as well. The most 
obvious thing to note is that the Log-Bicubic and Log-Spline 
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Figure 8. The four panels show the accuracy of shear re- 
covery with four different combinations of galaxy type (Regu- 
lar/Irregular) and PSF form (Wa/Wb) as denoted at the top of 
each panel. For the results in each panel, we use 10000 mock galax- 
ies, and scale them to four different sizes (galaxy radius equal to 
2.5, 5, 7.5, 10 times the PSF scale radius) to recovery the input 
shear values. In all panels, the PSF scale radius is half of the pixel 
size. The red and black data points are from the Log-Bicubic and 
Log-Spline interpolation methods respectively. The input shear 
values are shown by dashed lines. From the top to the bottom 
panel, the ratio of the mean flux density of the noise to that of 
the galaxy (In/Ig) is 0.1, 0.5, 0.6, and 0.2 respectively. 
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interpolation methods are accurate image reconstruction ap- 
proaches not only for weak lensing, but also for all kinds of 
other purposes. The way we remove the noise contamination 
from the shear signal can in principle also be considered in 
other shear measurements, in particular those that are based 
on measuring the multipole moments of the source images 
{e.g. , Kaiser et al. 1995 and its various extensions). 

So far, our discussion has neglected the "photon count- 
ing" shot noise, which is always present due to the finite 
telescope exposure time. Indeed, it becomes the dominant 
source of photon noise for ground-based weak lensing sur- 
vey because of the large sky background 2 . Unlike the as- 
tronomical photon noise that we have discussed, the pho- 
ton shot noise varies from pixel to pixel independent of the 
size of the PSF. Furthermore, the fluctuation amplitude of 
the shot noise is also dependent on the photon flux of the 
source. Therefore, it is hard to cleanly remove the contam- 
ination from the photon counting shot noise in our shear 
measurement. We do not intend to deal with this problem 
in this paper. Alternatively, we argue that the systematic 
shear measurement error due to the shot noise can be sup- 
pressed by simply increasing the telescope exposure time. 
More specifically, we argue that for any given tolerance level 
of the systematic error, there is a critical exposure time, be- 
yond which the contamination from the shot noise is ade- 
quately suppressed. To demonstrate this statement, in fig. 9, 
we plot the accuracy of shear recovery under three different 
conditions: high, medium, and low noise level, which cor- 
respond to noise-to-source mean surface brightness ratio of 
100, 10, and 1, or the signal-to-noise ratio of 0.01, 0.1, 1 
respectively. The x axis in each plot is the mean number 
of background noise photons recorded on each pixel, which 
is proportional to the exposure time. The accuracy of shear 
recovery is expressed in terms of the multiplicative and addi- 
tive bias parameters (commonly used by people in the weak 
lensing community), which are defined as: 

7measured = 7true(l + m) + C (11) 

For a perfect shear measurement, one should have m = c — 
0. To produce each (m, c) pair in fig. 9, we use five differ- 
ent input shear values, i.e. , (0.04, 0.04), (0.02, 0.02), (0, 
0), (-0.02, -0.02), (-0.04, -0.04), for (71, 72) (which are re- 
ported using red and blue colors respectively), and 10000 
mock irregular galaxies for each input (71, 72). Note that 
to see the dependence of the shear recovery accuracy on the 
exposure time more clearly, we repeatedly use the same set 
of galaxies for different exposure times. In all cases, we set 
the galaxy radius to be 7.5 times the PSF radius. The later 
is set to be equal to the pixel size. We use Log-Spline as the 
interpolation method. The astronomical photon noise is not 
included, i.e. , the photon counting noise is the only photon 
noise in this test. According to the figure, the systematic er- 
rors due to the photon counting noise is clearly suppressed 
when exposure time is beyond some threshold. Not surpris- 
ingly, the low noise case requires the shortest exposure time 



2 Note that the mean of the background photon does not affect 
our shear measurement, because the method of Z08 uses the spa- 
tial derivatives of the surface brightness field. It is the fluctuation 
of the photon numbers across the pixels that may bias our shear 
estimator. 



to achieve the same shear recovery accuracy level. A more 
comprehensive test with a much larger ensemble of galaxies 
will be shown in a separate paper. 

The other sources of systematic errors that we have ne- 
glected include the high order corrections (e.g. , 7 2 , k 2 , jk) 
to our master equation [eq.(5)], the spatial variations of the 
cosmic shear, etc.. These factors likely affect the measured 
shear values at percent levels on cosmic scales, which is im- 
portant in the era of precision cosmology. For clustering lens- 
ing, the high order shear terms are more important because 
the shear is of order ten percent. This subject will be studied 
in a companion paper. 

This paper is a natural continuation of Z08 on the 
methodology of cosmic shear measurement. In another pa- 
per, we will further test this method with the data from the 
Shear TEsting Program (Heymans et al. 2006; Massey et al. 
2007) and the GREAT08 program (Bridle et al. 2009), and 
also present results measured with real astronomical data. 
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APPENDIX - DEFINITIONS OF THE 
INTERPOLATION METHODS 

In this appendix, we give the mathematical definitions of the 
three classic 2D interpolation methods: Bilinear, Bicubic, 
Spline. For their logarithmic extensions [i.e. , Log-Bilinear, 
Log-Bicubic, Log-Spline), we only have one minor point to 
address at the end of this section. 

The Bilinear method is the simplest of the three. Let 
us write the coordinates of the grid points as (xi,y.j) (i,j — 
1,2,3...), and the signals as A(xi,yj). Suppose the point of 
our interest is (x,y), which satisfies Xi < x < x i+ i and 
Hi ^ U ^ tnc Bilinear method defines A(x,y) in the 

following way: 

A(x,y) = tuA(x l+1 ,y j+1 ) + (l-t)uA(xi,y J+1 ) (12) 
+ t(l - u)A(x i+1 , yj ) + (1 - t){l - u)A( Xi , yj ) 

where 



t = (x - Xi)/(x i+1 - Xi) 

u=(y- %')/(%■+! - Vj)- 



(13) 



The Bicubic method includes higher order terms of t 
and u to achieve smoothness of the interpolated function. It 
requires the user to specify not only the signal A(xi, yj), but 
also the spatial derivatives dA/dx, dA/dy, and d 2 A/dxdy 
at every grid point (xi,yj). Since the spatial derivatives of 
the signal are usually not known a priori, we estimate them 
using the finite-difference method: 



dx *' 3 Xi+i — Xi-i 

9A _ A(xj,y j+1 ) - A{x i ,y j - 1 ) 

oy yj+i - Vi-i 



(14) 



d A 
dxdy 



(xi,yj) = [A(x i+1 ,y j+1 ) + A(xi-i, yj -i) 

- A(x i+ i,yj-i) - A(x i - 1 ,y j+1 )] 
I [(x i+ i - Xi-i)(y j+ i - yj-i)] 



The interpolated function inside each grid square is written 
in the following polynomial form: 



A(x,y) = J2J2c mn t m - 1 u n - 1 



(15) 



The values of the sixteen parameters c mn are constrained 
using eq.(15) and the following three equations at the four 
corners of the grid square: 



^(*.») = EE(™-i)«w*^ a --"- 1 ' / ' 



in — '2 n=l 
4 4 



u 'Tx ^ 



dy 



( x j v) = E E^ n_ 1 ) c " 



it U 



dy 
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OA 

dxdy 



(x, y) = — - l)c m „i" 



m— 2 n — 2 



-2 u n-2 d« 

da; dj/ 



where f and w have been defined in eq.(13). 

The 2D Spline method simply refers to using the ID 
Spline interpolation along each dimension. The ID (cubic) 
Spline interpolation method works as follows: 

Given the values of a function f(x) at a set of points x% 
(i = 1...N), the form of the function in the interval between 
Xj and Xj+i is written as: 

f(x) = Hf{xi) + Kf(x J+1 ) + Uf"( Xj ) + Vf"(x J+1 ) (17) 



where 



H = 



(18) 



K = l-H 
V=^(K 3 -K)(x j+1 



and f"{xj) is the second derivative of the function / at Xj. 
As a consistency check, one can easily show that d 2 f /dx 2 = 
Hf"(xj) + K f" (xj+i). The value of the second derivatives 
are specified by requiring the first derivatives evaluated from 
the two sides of the grid point to be equal. Note that this 
requirement only provides iV— 2 equations, while there are N 
second derivatives in total. The rest of the constraint comes 
from the boundary conditions on /"(l) and f"(N). In this 
paper, we simply set /"(l) = f"(N) — 0, which yields the 
so-called natural cubic spline. 

Finally, we note that the "Log" based interpolation 
methods are all well defined except when the readouts of 
some pixels are zero. This is a very rare case in practice 
due to the presence of noise. However, this situation can in 
principle exist in simulations. To cure this problem, one can 
either change the zeros into tiny positive numbers, or simply 
avoid interpolating the regions with zeros. The second option 
says that if a grid square (regarding the Log-Bilinear and 
Log-Bicubic methods) or a unitary segment (regarding the 
Log-Spline method) contains any zero readouts in their four 
corners or two ends, the finer grid points within them are all 
set to have zero values. The rest of the grid squares/segments 
are interpolated independently as usual. Note that in the 
Log-Spline method, this means the Spline interpolations are 
carried out only in those nonzero segments that are isolated 
by the zeros. These two choices usually work similarly well. 
However, when there are extended regions of zero readouts, 
we find that the second option is better, because it avoids 
introducing artificial high order fluctuations in the zero re- 
gions by methods like Log-Bicubic or Log-Spline. 
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Figure 9. We show that in the shear measurement method of 
Z08, how the shear recovery accuracy is affected by the exposure 
time (i.e. , the mean number of background photons received per 
pixel). The results are expressed in terms of the multiplicative 
(m) and additive (c) bias parameters defined in eq. (11). Each 
data point is achieved using five different sets of (71, 72), i.e. , 
(0.04, 0.04), (0.02, 0.02), (0, 0), (-0.02, -0.02), (-0.04, -0.04). The 
red and blue data points are for 71 and 72 respectively. For each 
input shear, 10000 mock galaxies are used to recover the shear. 
The same galaxies are used repeatedly for different exposure times 
so that the dependence of the shear recovery accuracy on the 
exposure time can be more clearly seen. The exposure time is 
denoted by the mean number of background photons received 
per pixel (x-axis). The three panels from the top to the bottom 
correspond to the signal-to-noise ratios (i.e. , the ratio of the 
mean surface brightness of the source to that of the background 
noise) of 0.01, 0.1, 1 respectively, referring to the high, medium, 
and low noise cases as denoted at the top of the panels. 
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